#!/usr/bin/perl
package run_yn00;

require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(	
					YN00_KAKS
					YN00_KA__
					YN00_KS__
					YN00_A___
					YN00_S___
					prepare_calculate_yn00
					calculate_yn00
					cleanup_calculate_yn00
					);

use strict;
use warnings;
#use parse_bl2seq;

sub cleanup_calculate_yn00($)
{
	my ($dir) = @_;
	unlink(glob("$dir/yn00_tmp/*"));
	rmdir("$dir/yn00_tmp");
}

# for alignment results
use constant YN00_KAKS =>  0;
use constant YN00_KA__ =>  1;
use constant YN00_KS__ =>  2;
use constant YN00_A___ =>  3;
use constant YN00_S___ =>  4;

sub calculate_yn00($$$)
{
	my ($dir, $seq1, $seq2) = @_;
	die unless length $seq1 == length $seq2;
	open FH_SEQ, ">$dir/yn00_tmp/seq.tmp";
	print FH_SEQ "2 ".length($seq1), "\n";
	print FH_SEQ "a\n", $seq1, "\n";
	print FH_SEQ "b\n", $seq2, "\n";
	close FH_SEQ;
	unlink "$dir/yn00_tmp/yn";
	
	chdir "$dir/yn00_tmp";
	system ("/share/tools/paml/yn00 < yn00.ctl > /dev/null");
	chdir "..";
	unless (open FH_YN, "$dir/yn00_tmp/yn")
	{
		return;
	}
	while (<FH_YN>)
	{
		next unless /seq\.\s+seq./;
		<FH_YN>;
		$_ = <FH_YN>;
		my @values = split /[\s\+\-]+/;
		shift @values;
		return unless @values == 11;
		my @results = ($values[6],                 # kaks
					   $values[7],                 # ka
					   $values[9],                 # Ks
					   $values[3],                 # A
					   $values[2]                  # S
					   );
	#	print STDERR join "\n", @results, "\n";
		# ignore invalid values
		for (@results)
		{
			unless (defined $_ && /\d/ && $_ >= 0.0)
			{
				$_ = '\N' ;
			}
		}
		
		# ignore out of range
		for (YN00_KAKS, YN00_KA__, YN00_KS__)
		{
			$results[$_] = '\N'  unless ($results[$_] < 50.0);
		}
		return (@results);				
	}
	return;
}





1;

